3. Bölüm: R ile Çoklu Regresyon Analizi: (Wooldridge Örnekleri) Ders Notu

Önceki bölümlerde basit regresyon modeli ele alındı ancak basit regresyon açıklamak istediğimiz bağımlı değişkeni sadece bir bağımsız değişkenle açıklamaya çalıştığı için ceteris paribus (diğer herşey sabitken) varsayımını ihlal etmesi çok olasıydı. Üstelik bu bağımlı değişkeni açıklayan ve bağımlı değişkenin üzerinde de etkisi olan başka bir değişken olmaması olasılığı da çok düşüktü. Çoklu regresyon modeli tahmincimize başka değişkenler ekleyerek y nin üzerindeki varyasyonu daha fazla açıklama şansı doğuracaktır. Çoklu regresyon modelinin bir diğer avantajı bağımsız değişkenin derecelerini arttırarak parametrelerde doğrusal olmasına rağmen doğrusal olmayan fonksiyonlar elde edilebilir ve daha iyi bir açıklayıcı elde edilebilir. Çoklu regresyon analizinde bir çok fonksiyonel form kullanılabilir ve bu bize daha fazla esneklik sağlar.

İki bağımsız değişkenli model

Bu bölümün Wooldridge kitabınızda bulunan ilk örnek maaş (wage) bağımlı değişkenin iki bağımsız değişkenle açıklamaya çalıştığımız model. Bu modeli tekrar hatırlayalım.

\[ wage= \beta_0 + \beta_1 educ + \beta_2 exper + u \]

Burada $exper$ iş hayatında ki deneyim yılı, $educ$ eğitim yılını gösteren bağımsız değişkenler, ve $u$ ise gözlemleyemediğimiz diğer faktörler olarak özetlenebilir. Burada basit regresyon modelinde de anlamaya çalıştığımız gibi asıl ilgilendiğimiz parametre $\beta_1$ yani maaşı etkileyen diğer bütün faktörleri sabit tuttuğumuzda eğitimin gelir üzerinde ne kadar bir etkisi olduğu.

Basit regresyon bölümünde kurduğumuz modeli dikkate alalım

\[ wage= \beta_0 + \beta_1 educ + u \] Burada $\beta_2 exper$ gözlemleyemediğimiz diğer faktörler $u$’nun içindeydi. Çoklu regresyon modelimizde $\beta_2 exper$ terimini $u$’nun içinden çıkarıp modelimize ekleyince en azından deneyimin ($exper$), maaşlar ($wage$) üzerindeki etkisini $\beta_2$ katsayısıyla ölçebileceğiz.

Ancak hala yapmamız gereken önemli varsayımlar mevcut. Bu ikili değişken modelinin bize kattığı artık $\beta_1$ katsayısını en azından deneyimi ($exper$) sabit tutarak ölçebileceğimizi biliyoruz. Basit regresyon modelinde, yani deneyim ($exper$), $u$ teriminin içindeyken, deneyimin, eğitimle bir korelasyonu olmadığını varsaymak zorundaydık. Burada şunu unutmamak gerekir, iki açıklayıcı değişkenli modelde de hala $u$ teriminin içinde bulunan gözlemleyemediğimiz değişkenler, eğitim ($educ$) ve deneyim ($exper$) değişkenleriyle ilgisiz olmalıdır. Bunun sağlanması çok zor olan bir varsayım olduğunu aklınızdan çıkarmayın.

Anlaşıldığı gibi, bu bölüm atlanan (gözlemlenemeyen) değişken sapmasını (“omitted variable bias”) da belli bir ölçüde açıklayacaktır. Çoklu regresyon modelini kullanmamızın bir diğer nedeni de örneğimizden anlaşılacağı üzere gözlemlenemeyen (atlanan) değişken sapmasını ortadan kaldırabilmektir.

Çoklu regresyonun ana fikirlerinden biri, bu ihmal edilen değişkenler hakkında verimiz varsa, onları ek regresörler olarak modele dahil edebilmek ve böylece diğer değişkenleri (deneyim ($exper$) gibi) sabit tutarken bir diğer regresörün (eğitim ($educ$) gibi) bağımsız değişken (maaşlar ($wage$) gibi) üzerindeki nedensel etkisini tahmin edebilmektir.

Alternatif olarak, eğer biri nedensel çıkarımlar değil de tahminle ilgileniyorsa, çoklu regresyon modeli, tek bir regresör kullanılarak yapılan tahminleri geliştirmek için birden fazla değişkenin regresör olarak kullanılmasını mümkün kılar. Bu konuya ilerideki kesimlerde daha detaylı olarak geri döneceğiz.

Kitabınızda sunulan bir diğer iki bağımsız değişkenli regresyon örneği, öğrenci başına yapılan harcamanın ($expend$), lise öğrencilerinin ortalama test sonuçları ($avgscore$) üzerindeki etkisini inceliyor. Örnek olarak şu model kullanılıyor.

\[ avgscore= \beta_0 + \beta_1 expend + \beta_2 avginc + u \] $avginc$ ailenin ortalama geliri ve $u$ diğer gözlemlenemeyen faktörleri temsil eder. Bulmaya çalıştığımız etki, öğrenci başına yapılan harcamanın ($expend$), lise öğrencilerinin ortalama test sonuçları ($avgscore$) üzerindeki etkisini gösteren $\beta_1$ dir. Ancak ailenin ortalama geliri, harcamalar üzerinde bir etkiye sahip olduğundan, $avginc$ değişkenini regresyona dahil etmemek, en küçük kareler metoduyla bulduğumuz basit regresyon katsayısını ($\beta_1$) sapmalı hale getirecektir.

İki bağımsız değişkenli çoklu regresyon modelini genel bir şekilde yazarsak,

\[ y= \beta_0 + \beta_1 x_1+ \beta_2 x_2 + u \]

Burada,

$\beta_0$ kesen, $\beta_1$ diğer faktörleri sabit tutarak y’deki değişimi x1’e göre ölçer $\beta_2$ diğer faktörleri sabit tutarak y’deki değişimi x2’ye göre ölçer

İki açıklayıcı değişkenli regresyon modelinin kitabınızdaki bir diğer örneği, doğrusal regresyon modelinin fonksiyonel olarak yazılabildiğini göstermek için verilmiştir. Ailenin tüketimini, ailenin gelirinin ikinci dereceden (quadratic) bir fonksiyonu olarak gösterirsek

\[ cons= \beta_0 + \beta_1 inc+ \beta_2 inc^2 + u \]

$cons$ ailenin tüketimi, $inc$ ailenin geliri, $u$ tüketimi etkileyen diğer faktörler olarak özetlenebilir.

k bağımsız değişkenli model

Çoklu regresyon modeli olarak adlandırılabilir.

\[ y= \beta_0 + \beta_1 x_1 + \beta_2 x_2 + \beta_3 x_3 + ... + \beta_k x_k + u \]

$u$ terimi hata terimi veya gözlenemeyen diğer faktörler olarak adlandırılabilir.

Çoklu regresyon modelinin temel varsayımını hatırlayalım.

\[ E(u | x_1, x_2, x_3,...,x_k)=0 \] Gözleyemediğimiz diğer faktörlerin, bağımsız değişkenlerle ilişkisiz olması gerektiğini gösteren varsayımdır.

Örnek3.1 Kolej Gpa’sının Belirleyicileri

ödeviniz olarak da verdiğim kolej Gpa’sının belirleyicileri örneğinin datasını indirelim ve içinde bulunan değişkenlere göz atalım.

library(wooldridge)
data(gpa1)

library(rmarkdown)
paged_table(gpa1)

Örnek için kullanacağımız değişkenlerin tanımlarını şöyle yapabiliriz.

colGPA: MSU GPA, Michigan State Üniversitesi not ortalaması

hsGPA: high school GPA, lise not ortalaması

ACT: ‘achievement’ score, ACT, bir öğrencinin okulda ne öğrendiğini ölçen bir başarı testidir. SAT daha çok bir yetenek testidir, muhakeme ve sözel yetenekleri test eder. Üniversiteler SAT veya ACT sonuçları ile öğrenci kabul eder.

Kitabınızda regresyon sonuçları direk verilen örneğin sonuçlarını bulmak isterseniz şu komutları kullanabilirsiniz.

coklureg <- lm(colGPA~ hsGPA+ACT,data = gpa1)
summary(coklureg)

## 
## Call:
## lm(formula = colGPA ~ hsGPA + ACT, data = gpa1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85442 -0.24666 -0.02614  0.28127  0.85357 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.286328   0.340822   3.774 0.000238 ***
## hsGPA       0.453456   0.095813   4.733 5.42e-06 ***
## ACT         0.009426   0.010777   0.875 0.383297    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3403 on 138 degrees of freedom
## Multiple R-squared:  0.1764, Adjusted R-squared:  0.1645 
## F-statistic: 14.78 on 2 and 138 DF,  p-value: 1.526e-06

\[ colGPA= 1.29 + 0.453 hsGPA+ 0.0094 ACT + u \]

Peki bu modeli nasıl yorumlayacağız? İlk olarak kesişim katsayısı, hsGPA ve ACT’nin her ikisi de sıfır olursa, öngörülen üniversite GPA’sini göstermektedir ve örneğimizde 1.29 olarak hesaplanmıştır. Üniversiteye giden hiç kimsenin lise not ortalaması sıfır veya başarı testi sıfır olmadığı için, bu denklemdeki kesişim kendi başına anlamlı değildir. Eğim katsayılarına bakınca ikisinin de bekleneceği gibi pozitif olduğunu görüyoruz. Eğer ACT puanını sabit tutarsak, hsGPA’deki 1 birimlik artış colGPA üzerinde 0.453 puanlık bir artış yaratacaktır. Başka bir deyişle, eğer ACT puanları aynı olan iki öğrenci varsa ve bir öğrencinin lise not ortalaması diğer öğrenciden bir puan yüksekse, bu öğrencinin üniversite not ortalamasının yarım puan daha fazla olmasını tahmin edebiliriz.

Üniversite not ortalamalarının 4 üzerinden, ACT sonuçlarının 36 puan üzerinden hesaplandığı düşünülürse, ACT’nin eğimi 0.0094’ün üniversite not ortalaması üzerinde fazla bir etkisi olmadığı düşünülebilir. Hata paylarını göz önüne aldığımızda ACT’nin sadece küçük değil aynı zamanda istatiksel olarak da anlamlı olmadığını göreceğiz.

Basit regresyon örneğimizi hatırlarsak

basitreg <- lm(colGPA~ACT,data = gpa1)
summary(basitreg)

## 
## Call:
## lm(formula = colGPA ~ ACT, data = gpa1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.85251 -0.25251 -0.04426  0.26400  0.89336 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.40298    0.26420   9.095  8.8e-16 ***
## ACT          0.02706    0.01086   2.491   0.0139 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3656 on 139 degrees of freedom
## Multiple R-squared:  0.04275,    Adjusted R-squared:  0.03586 
## F-statistic: 6.207 on 1 and 139 DF,  p-value: 0.0139

ACT’nin katsayısı çoklu regresyona göre bir hayli büyük kalıyor (0.02706/0.009426=2.87 kat daha fazla). ACT değerinin çoklu regresyonda lise not ortalaması eklenince azaldığını gözlemliyoruz. Çünkü lise not ortalaması ve ACT arasında yüksek bir korelasyon var. Bu durum bize aslında basit regresyon kullanıldığında temel varsayımı ihlal ettiğimize dair güçlü bir ipucu veriyor.

Bu örnekte son olarak basit ve iki bağımsız değişkenli regresyon grafiklerini gösterelim. Basit regresyon örneğimiz için, dağılım grafiğine bir regresyon çizgisi eklememizin yeterli olduğunu hatırlayalım.

plot(gpa1$ACT,gpa1$colGPA)
abline(basitreg)

Bu grafik sadece iki boyutlu olduğundan, iki boyutlu bilgisayar ekranlarımızda göstermek daha kolaydır. Y eksenine üniversite not ortalaması, X eksenine ACT puanlarını koyuyoruz. Son olarak bulduğumuz eğim ve kesen katsayılarını grafiğimizde temsil eden doğrumuzu çiziyoruz.

Ancak çok değişkenli regresyon grafiklerini, iki boyutlu bilgisayar ekranlarda göstermek biraz daha zor olabilir. En azından örneğimizdeki gibi üç boyutlu bir grafiğe ihtiyacımız var. Önce bütün değişkenler için üç boyutlu bir dağılım grafiği çizelim daha sonra doğru yerine bir düzlem eklememiz gerekecek.

library(plotly)
library(reshape2)
plot_ly(data = gpa1, z = ~colGPA, x = ~hsGPA, y = ~ACT, opacity = 0.5) %>%
  add_markers()

Üç boyutlu grafiğimiz olduğu için artık üzerine gidip sağa sola çevirebilirsiniz. Tahmin edebileceğiniz gibi üç boyutlu grafikleri algılamamız daha büyük boyutlu grafiklerden daha kolay. Şimdi bu noktalarda oluşacak hata payını minimize edecek bir düzleme ihtiyacımız var. Bunun için çoklu regresyonumuzdan yardım alacağız.

cf.mod <- coef(coklureg)
hsGPA <- seq(min(gpa1$hsGPA),max(gpa1$hsGPA),length.out=50)
ACT <- seq(min(gpa1$ACT),max(gpa1$ACT),length.out=50)
colGPA <- t(outer(hsGPA, ACT, function(x,y) cf.mod[1]+cf.mod[2]*x+cf.mod[3]*y))
p <- plot_ly(x=~hsGPA, y=~ACT, z=~colGPA,type="surface") %>%
  add_trace(data=gpa1, x=gpa1$hsGPA, y=gpa1$ACT, z=gpa1$colGPA, mode="markers", type="scatter3d", opacity=0.5)
p

Sizin şimdilik bu grafikleri çizerken kullandığımız kodları anlamanıza gerek yok. Sadece grafiklerin şekline odaklanın. Bu grafik bize yapmış olduğumuz minimizasyon işleminin neden basit regresyona bir değişken daha eklediğimizde farklılaştığı ve ek değişkenin toplamda neden daha çok varyasyon açıklamaya yardım ettiği hakkında daha fazla bir içgüdü sunabilir.

Regresyon Tahmin Değerleri ve Hata Payı (OLS Fitted Values and Residuals)

Kitabınızdaki bu bölümde çoklu regresyonda tahmin değerleri ve hata payı üzerinde duruluyor. Tahmin değerlerini şapkayla (hat) gösteriyorduk.

\[ \hat{y_i}= \hat{\beta_0} + \hat{\beta_1} x_{i1} + \hat{\beta_2} x_{i2} + \hat{\beta_3} x_{i3} + ... + \hat{\beta_k} x_{ik} \] Gerçek popülasyon modelimiz ise bu şekilde gösteriliyordu.

\[ y_i= \beta_0 + \beta_1 x_{i1} + \beta_2 x_{i2} + \beta_3 x_{i3} + ... + \beta_k x_{ik} + u_i \] Normalde, bir $i$ gözlemi için gerçek değer $y_i$ ve tahmin edilen $\hat{y_i}$ birbirine eşit olmaz. Biz ancak bu hata payını minimize eden en küçük kareler yöntemini kullanarak doğrusal regresyon tahminlerini kullanabiliriz. Bu durumda her bir gözlem için hata payını şu şekilde tahmin edebiliriz.

\[ y_i-\hat{y_i}=\hat{u_i} \]

Çoklu regresyon tahmin değerleri ve hata payının özellikleri.

Hata payının örnek ortalaması sıfırdır, bu yüzden $\bar{\hat{y}}=\bar{y}$ olur.
Hata payının herbir bağımsız değişkenle örnek kovaryası sıfırdır. Bu yüzden çoklu regresyon tahmin değerleri ve hata payı arasında örnek kovaryans da sıfırdır.
Bütün bağımsız değişkenlerin ve bağımlı değişkenin ortalamaları regresyon çizgisi üzerindedir.

\[ \bar{y}= \hat{\beta_0} + \hat{\beta_1} \bar{x_1} + \hat{\beta_2} \bar{x_2} + \hat{\beta_3} \bar{x_3} + ... + \hat{\beta_k} \bar{x_k} \]

Çoklu Regresyon Modeli (Multiple Linear Regression (MLR) Model) Varsayımları

Birinci Varsayım (MLR.1): Parametrelerde doğrusaldır.

\[ y= \beta_0 + \beta_1 x_{1} + \beta_2 x_{2} + \beta_3 x_{3} + ... + \beta_k x_{k} + u \] $\beta$ terimleri doğrusaldır ancak $x$ terimlerinin doğrusal olmasına gerek yoktur. Çoğu zaman $x$ terimlerinin logaritması alınabilir veya $x^2$ gibi fonksiyonal formları kullanılabilir.

İkinci Varsayım (MLR.2): Rastgele Örnekleme.

Rastgele n örnek gözlemimiz var.

Üçüncü Varsayım (MLR.3): Mükemmel bir eş doğrusallık yoktur.

Örneklemde (ve dolayısıyla popülasyonda), bağımsız değişkenlerin hiçbiri sabit değildir ve bağımsız değişkenler arasında kesin doğrusal ilişkiler yoktur.

Dördüncü Varsayım (MLR.4): Sıfır koşullu ortalama.

Modelde bulunan bağımısız değişkenlerin herhangi bir değeri için hata payı $u$’nun koşullu beklenen değeri sıfırdır. $E(u | x_1, x_2, ..., x_k) = 0$

Gözlemlenemeyen Değişken Sapması (Omitted Variable Bias)

Gerçek popülasyon modelinde bulunan bir değişkeni modelimize eklememek gözlenemeyen değişken sapmasına yol açabilir. Bu sapmayı kitabınızda verilen örnek üzerinden inceleyeceğiz.

Varsayalım ki maaşları açıklayabilen sadece iki değişkene ihtiyacımız var. Bunlar eğitim ve yetenek. O zaman gerçek popülasyon modelimiz şu şekilde olur.

\[ maaş = \beta_0 + \beta_1 eğitim + \beta_2 yetenek + u \]

Bu modelin bütün çoklu regresyon varsayımlarını (MLR1-MLR4) sağladığını varsayalım. Şimdi yeteneği gözlemleyemediğinizi varsayalım ve sadece, eğitim ve maaş verilerine sahip olduğunuzu düşünün. Sadece basit regresyon modelini kullanabilirsiniz ve $\beta$ katsayılarını tahmin edebilirsiniz. Burda gerçek modeli kullanamadığınız için bu regresyon modelinden elde ettiğimiz tahminleri tilda $\sim$ ile gösterelim.

\[ maaş = \beta_0 + \beta_1 eğitim + v \]

\[ \tilde{maaş} = \tilde{\beta_0} + \tilde{\beta_1} eğitim \] $v$ hata payı, yetenek gözlenemediği için şu şekilde yazılabilir: $v=\beta_2 yetenek + u$

Eğer $yetenek$ gözlemlenebilseydi şu şekilde tahmin edebilirdik.

\[ \hat{maaş} = \hat{\beta_0} + \hat{\beta_1} eğitim + \hat{\beta_2} yetenek \] $\hat{\beta_1}$ tahmincisi yetenek değişkenini sabit tutarak bir birim eğitim değişimin maaşı kaç birim değiştireceğini göstermektedir. $\tilde{\beta_1}$ ise bir birim eğitim değişimin maaşı kaç birim değiştireceğini gösterir ancak yetenek değişkenini sabit tutmaz.

Burdan $\tilde{\beta_1}$ sapması şu şekide gösterilebilir

\[ \tilde{\beta_1}= \hat{\beta_1} + \hat{\beta_2} \tilde{\delta_1} \] $\tilde{\delta_1}$ yeteneğin eğitim üzerindeki basit regresyonunun eğim parametresidir.

\[ \tilde{eğitim} = \tilde{\delta_0} + \tilde{\delta_1} yetenek \] Eğitim yetenek arttığı için artacağından, yeteneği modele katmamak $\tilde{\beta_1}$ tahmininin o kadar artmasına yol açar.

Şimdi bu hesaplamamızı bir örnek üzerinden gösterelim. Bu örneği hatırlarsak

##install.packages("stargazer")
library(stargazer)

## 
## Please cite as:

##  Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.

##  R package version 5.2.2. https://CRAN.R-project.org/package=stargazer

stargazer(list(coklureg,basitreg),type = "text")

## 
## =================================================================
##                                  Dependent variable:             
##                     ---------------------------------------------
##                                        colGPA                    
##                               (1)                    (2)         
## -----------------------------------------------------------------
## hsGPA                      0.453***                              
##                             (0.096)                              
##                                                                  
## ACT                          0.009                 0.027**       
##                             (0.011)                (0.011)       
##                                                                  
## Constant                   1.286***               2.403***       
##                             (0.341)                (0.264)       
##                                                                  
## -----------------------------------------------------------------
## Observations                  141                    141         
## R2                           0.176                  0.043        
## Adjusted R2                  0.164                  0.036        
## Residual Std. Error    0.340 (df = 138)       0.366 (df = 139)   
## F Statistic         14.781*** (df = 2; 138) 6.207** (df = 1; 139)
## =================================================================
## Note:                                 *p<0.1; **p<0.05; ***p<0.01

Çoklu regresyon modeli şu şekilde yazılabilir

\[ colGPA = \beta_0 + \beta_1 ACT + \beta_2 hsGPA + u \] Çoklu regresyonun sonuçları şu şekilde rapor edilebilir

\[ colGPA = 1.286 + 0.009 ACT + 0.453 hsGPA \] Basit regresyon modeli

\[ colGPA = \beta_0 + \beta_1 ACT + v \] Basit regresyon modeli sonuçları

\[ colGPA = 2.403 + 0.027 ACT \] Bu sonuçlar sapmalı sonuçlar olduğu için $\sim$ işaretiyle gösterelim

\[ \tilde{colGPA} = \tilde{\beta_0} + \tilde{\beta_1} ACT \]

Çoklu regresyon modelimiz gerçeğe daha yakın bir model olduğundan tahminleri şapkayla gösteriyoruz.

\[ \hat{colGPA} = \hat{\beta_0} + \hat{\beta_1} ACT + \hat{\beta_2} hsGPA \] Sapmalı tahminleri bu formülle gösteriyorduk

\[ \tilde{\beta_1}= \hat{\beta_1} + \hat{\beta_2} \tilde{\delta_1} \] Bütün bildiğimiz değerleri yerine koyarsak

\[ 0.027= 0.009 + 0.453 \tilde{\delta_1} \] Önce işlemle sonra basit regresyonla sapmayı $\tilde{\delta_1}$ bulalım.

\[ 0.027= 0.009 + 0.453 \tilde{\delta_1}\\ \tilde{\delta_1}= \frac{(0.027-0.009)}{0.453}=0.0397 \] $\tilde{\delta_1}$ iki bağımsız değişkenin birbirleri arasındaki basit regresyon ilişkisinden de bulunabilir

\[ \tilde{ACT} = \tilde{\delta_0} + \tilde{\delta_1} hsGPA \]

sapmareg<-lm(hsGPA~ACT,data = gpa1)
stargazer(list(sapmareg),type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                hsGPA           
## -----------------------------------------------
## ACT                          0.039***          
##                               (0.009)          
##                                                
## Constant                     2.463***          
##                               (0.218)          
##                                                
## -----------------------------------------------
## Observations                    141            
## R2                             0.120           
## Adjusted R2                    0.113           
## Residual Std. Error      0.301 (df = 139)      
## F Statistic           18.879*** (df = 1; 139)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Regresyon sonucuyla bulduğumuz ve hesapladığımız $\tilde{\delta_1}$ aynı olduğuna dikkat edin. Basit regresyon yaptığımızda gözlemleyemediğimiz değişken yüzünden olan sapma $\hat{\beta_2} \tilde{\delta_1}=0.453 \cdot 0.0397=0.017$ olur. Bu fark $\tilde{\beta_1}-\hat{\beta_1}$ arasındaki farka eşittir.

Çoklu regresyon (En küçük kareler (OLS)) tahmincilerinin varyansı

Homoskedastisite (Eş varyanslık), Çoklu eşdoğrusallık (Multicollinearity) ve Varyans enflasyon faktörü (Variance inflation factor (VIF))

Çoklu regresyonun 5. varsayımını hatırlayalım.

Beşinci Varsayım (MLR.5): Hata terimi $u$ açıklayıcı değişkenlerin her bir değeri için aynı varyansa sahiptir.

MLR.1 ila MLR.5 arasındaki varsayımlar topluca Gauss-Markov varsayımları olarak bilinir (kesitsel regresyon için).

Bir eğimin varyansını hesaplamak için şu formül kullanılabilir.

\[ var(\hat{\beta_j})=\frac{\sigma^2}{SST_j(1-R_j^2)} \] $SST_j$’yi hatırlarsak

\[ SST_j = \sum{(x_{ji}-\bar{x_j})^2}=(n-1)(var(x_j)) \] Bu durumda $var(\hat{\beta_j})$’yi üç parça halinde yazabiliriz.

\[ var(\hat{\beta_j})=\frac{1}{n}\cdot\frac{\sigma^2}{var(x_j)}\cdot\frac{1}{1-R_j^2} \]

Bu formülle birlikte şu çıkarımları yapabiliriz.

$\frac{1}{n}$ bize n yani gözlem sayısı arttıkça eğimin varyansının azalacağını söyler.
$\sigma^2$ bize hata teriminin çok fazla değişiklik gösterdiğinde, $\beta$ varyansının artacağını söyler.
$\frac{1}{var(x_j)}$ bize $x_j$’nin varyansı arttıkça $\beta$ varyansının azalacağını söyler. Bunun nedeni $x_j$’in bu sayede konuyla ilgili daha ilişkili bilgiler sağlayabilecek olmasıdır.
$\frac{1}{1-R_j^2}$ varyans enflasyon faktörü (Variance inflation factor (VIF)) olarak adlandırılır ve çoklu eşdoğrusallık (Multicollinearity) ölçer. Eğer $x_j$ bağımlı değişkeni diğer regresörlerle çok ilişkiliyse (korelasyon yüksekse) bu $R_j^2$’yi arttırır. $R_j^2$ bire yaklaştıkça $\frac{1}{1-R_j^2}$ paydası sıfıra yaklaşır ve $\frac{1}{1-R_j^2}$ sonsuza yaklaşır. Dolayısıyla $\beta_j$’in varyası çok fazla artar.

Örneğimizi kullanarak bu parçaları çıkaralım.

standarthata<-summary(coklureg)$sigma
standarthata

## [1] 0.3403158

Standart hatayi hesaplamak için sadece sigma komutunu kullanıyoruz. Standart hatayi colGPA örneği için 0.34 olarak hesapladık. $\sigma^2$ yani standart hatanın karesi bize regresyonumuzun hata payının varyansını verecektir.

sigmakare<-standarthata^2
sigmakare

## [1] 0.1158148

Regresörlerimiz arasındaki $R_j^2$ hesaplamak için daha önce hesapladığımız sapmareg’i kullanabiliriz. Bu $R_j^2$ bizim için lise not ortalamasıyla ACT sonuçları arasındaki ilişkiyi gösterecektir. Eğer $R_j^2$ bir’e yakınsa çoklu eş doğrusallık muhtemel olabilir.

Rjkare<-summary(sapmareg)$r.squared
Rjkare

## [1] 0.1195815

İki bağımsız değişkenimizin basit regresyonundan elde ettiğimiz $R_j^2$ yaklaşık olarak yüzde 12’dir.

$\frac{1}{1-R_j^2}$ varyans enflasyon faktörü (Variance inflation factor (VIF))’i hesaplayalım.

VIF<-1/(1-Rjkare)
VIF

## [1] 1.135823

son olarak n’i bulmalıyız.

n<-nobs(coklureg)
n

## [1] 141

son olarak $var(x_j)$’yi bulmamız lazım.

varxj<-var(gpa1$hsGPA)*(n-1)/n
varxj

## [1] 0.1016267

Şimdi formülümüzü kullanarak $var(\hat{\beta_j})$’yi hesaplayabiliriz.

\[ var(\hat{\beta_j})=\frac{1}{n}\cdot\frac{\sigma^2}{var(x_j)}\cdot\frac{1}{1-R_j^2} \]

varbetaj<-(1/n)*(sigmakare/varxj)*VIF
varbetaj

## [1] 0.009180115

$\beta_j$’nin standart hatası $var(\hat{\beta_j})$’nin karaköküdür.

sdbetaj<-varbetaj^0.5
sdbetaj

## [1] 0.09581292

Tablodan hsGPA’in standart hatasına bakarsanız bulduğumuz rakamla aynı olduğunu göreceksiniz.

VIF’i bulmanın kolay bir yolu car paketini yüklemektir.

require(car)

## Loading required package: car

## Loading required package: carData

vif(coklureg)

##    hsGPA      ACT 
## 1.135823 1.135823

require komutu paket yükleme komutu (install.packages) ve library komutunun ortak komutu olarak düşünülebilir. Eğer bir paketi daha önce yükleyip yüklemediğinizi hatırlamıyorsanız. require paket yüklenmemişse yükleyecek ve library() yapacaktır. Eğer daha önceden yüklüyse sadece library() yapacaktır. vif() komutu VIF’i regresyonunuz için otomatik olarak hesaplar. 1.1358’in bizim daha önce hesapladığımız VIF değeriyle aynı olduğuna dikkat edin.

Ödev

Wooldridge Bölüm Soruları

4137 üniversite öğrencisi üzerinde GPA2’deki verileri kullanarak, en küçük kareler (OLS) tarafından aşağıdaki denklem tahmin edilmiştir:

\[ üniversite.not.ortalaması=1.392 - 0.0135 lmsyd + 0.00148sat \]

lmsyd, öğrencilerin lise mezuniyet sınıflarındaki yüzdelik dilimi göstermektedir. Örneğin, lmsyd=5, sınıfın ilk %5’inde olduğu anlamına gelir. sat, öğrenci başarı testindeki birleştirilmiş matematik ve sözel puanlardır.

lmsyd katsayısının negatif olması neden mantıklıdır?
lmsyd = 20 ve sat = 1.050 olduğunda tahmini üniversite not ortalaması nedir?
Diyelim ki iki lise mezunu, A ve B, liseden aynı yüzdelik dilimde mezun oldular. Ancak Öğrenci A’nın SAT puanı 140 puan daha yüksek. (aynı zamanda yaklaşık bir standart sapmaya sahip). Bu iki öğrenci için üniversite not ortalamasında tahmini fark nedir? Bu fark büyük bir fark mı?
lmsyd sabit tutulduğunda, SAT puanlarındaki hangi fark, tahmini üniversite not ortalaması farkına 0,50 veya bir notun yarısı kadar farka yol açar? Cevabınıza yorum yapın.

Aşağıdaki model, Biddle ve Hamermesh (1990) tarafından uyumak ve çalışmak için harcanan zaman arasındaki dengeyi incelemek ve uykuyu etkileyen diğer faktörleri incelemek için kullanılan çoklu regresyon modelinin basitleştirilmiş bir versiyonudur:

\[ uyku=\beta_0 + \beta_1 toplam.iş+ \beta_2 eğitim + \beta3 yaş + u \]

uyku ve toplam iş dakikayla, eğitim ve yaş yılla hesaplanmıştır.

Yetişkinler iş için uykuyu takas ediyorsa, $\beta_1$’in işareti nedir?
$\beta_2$ ve $\beta_3$’ün hangi işaretleri olacağını düşünüyorsunuz?
SLEEP75’teki verileri kullanarak, tahmin edilen denklem şu şekildedir:

\[ uyku=3,638.25 - 0.148 toplam.iş - 11.13 eğitim + 2.20 yaş \] $n=706$ ve $R^2=0.113$

Birisi haftada beş saat daha fazla çalışırsa, uykunun kaç dakika düşeceği tahmin edilir? bu büyük bir takas mı? d.Eğitim üzerindeki tahmini katsayının üzerindeki işaretini ve büyüklüğünü tartışın. e.Uykudaki çeşitliliğin çoğunu toplam iş, eğitim ve yaşın açıkladığını söyleyebilir misiniz? Uyuyarak geçirilen süreyi başka hangi faktörler etkileyebilir? Bunların toplam iş ile ilişkili olması muhtemel mi?

Üniversite not ortalamasını çeşitli etkinliklerde harcanan zamanla ilişkilendiren bir çalışma yapmak istiyorsunuz, birkaç öğrenciye bir anket dağıttınız. Öğrencilere her hafta dört aktivitede kaç saat geçirdikleri sordunuz: ders çalışmak, uyumak, bir işte çalışmak ve boş zaman. Herhangi bir aktivite dört kategoriden birine konur, böylece her öğrenci için dört aktivitedeki saatlerin toplamı 168 olmalıdır.

\[ not.ortalaması=\beta_0 + \beta_1 ders.çalışma+ \beta_2 uyumak + \beta3 iş.çalışma + \beta4 boş.zaman + u \]

$\beta_1$’i yorumladığınızda çalışmayı değiştirirken uykuyu, işi ve boş zamanları sabit tutmak mantıklı mı?
Bu modelin neden MLR.3 Varsayımını ihlal ettiğini açıklayın.
Modeli, parametrelerinin faydalı bir yoruma sahip olması için nasıl yeniden formüle edebilirsiniz ve varsayım MLR.3’ü ihlal edilmez?

Cevaplar

lmsyd, ne kadar küçükse öğrencinin lisedeki durumu o kadar düşük olacak şekilde tanımlanır. Diğer her şey eşittir, öğrencinin lisedeki durumu ne kadar kötüyse, beklenen üniversite not ortalaması o kadar düşük olur.
Sadece bu değerleri denkleme yerleştirin

\[ üniversite.not.ortalaması=1.392 - 0.0135 \cdot 20 + 0.00148 \cdot 1050 = 2.676 \]

A ve B arasındaki fark, sat katsayısının 140 katıdır, çünkü lmsyd her iki öğrenci için de aynıdır. Dolayısıyla A’nın $0.00148 \cdot (140)= 0.207$ daha yüksek bir puana sahip olduğu tahmin edilmektedir.
lmsy sabit olduğunda, $\sum(\Delta universite.not.ortalaması) = 0.00148 \cdot \Delta sat$. Bu durumda $0.5 = 0.00148 \cdot \Delta sat$ veya $sat= $ olacak şekilde bulmak istiyoruz. $sat= 338 $

Diğer herşey eşit olduğunda, yetişkinler uykuyu iş için tercih ediyorsa, daha fazla iş daha az uyku anlamına gelir, bu yüzden $\beta_1<0$ olur.
$\beta_2$ ve $\beta_3$ işaretleri kişiden kişiye göre değişir. Vereceğiniz örneklere göre bu katsayıların işaretleri eksi veya artı olabilir. Bu yüzden bu çalışma için herhangi bir beklenti içinde olamayız. Kimi yaş ilerledikçe uykunun azaldığını söyler, kimi arttığını ve bunu gerekçeleriyle örneklendirebilir. Aynı durum eğitim ile ilgili olarak da tartışılabilir.
Uyku ve toplam iş dakika ile ölçüldüğünden 5 saati dakikaya çevirmeniz gerekir ($5 \cdot 60= 300$). 300 dakikayı formülde yerine koyarsak $0.148 \cdot 300= 44.4$ dakika. Uyku bir hafta içinde 44.4 dakika düşer bunun çok fazla bir düşüş olduğunu söyleyemeyiz. d.Daha fazla eğitim, daha az tahmini uyku süresi anlamına gelir, ancak etkisi oldukça küçüktür. Üniversite ile lise arasındaki farkın dört yıl olduğunu varsayarsak, modelimiz diğer değişkenler sabit olduğunda üniversite mezunu olan kişinin lise mezunu olan bir kişiye göre haftada yaklaşık 45 dakika (11.13 ) daha az uyuduğunu tahmin ediyor.
Şaşırtıcı olmayan bir şekilde, üç açıklayıcı değişken, katılımcıların yalnızca yaklaşık %11.3’ünü ($R^2$) açıklamaktadır. Hata terimindeki önemli bir faktör genel sağlıktır. Bir diğeri ise medeni durum ve kişinin çocuğu olup olmadığıdır. Sağlık, medeni durum ve çocukların sayısı ve yaşları genellikle toplam iş süresiyle ile ilişkilidir.

Hayır. Tanım olarak, çalışma + uyku + iş + boş zaman = 168. Bu nedenle, çalışmayı değiştirirsek, toplamın hala 168 olması için diğer kategorilerden en az birini değiştirmeliyiz.
1. kısmından, mesela iş değişkenini diğer bağımsız değişkenlerin mükemmel bir lineer fonksiyonu olarak yazabiliriz: $iş = 168 - uyku - ders - boş zaman$ Bu her gözlem için geçerlidir, dolayısıyla MLR.3 ihlal edilmiştir.
Bağımsız değişkenlerden birini modelden çıkarın, mesela boş zaman değişkenini çıkarın:

\[ not.ortalaması=\beta_0 + \beta_1 ders.çalışma+ \beta_2 uyumak + \beta3 iş.çalışma + u \]

Bu durumda, $\beta_1$ ders çalışma bir saat arttığında, uyku, iş ve $u$ sabit tutulduğu zaman not ortalamasındaki değişiklik olarak yorumlanır. Uykuyu ve çalışmayı sabit tutuyoruz, ancak ders çalışmayı bir saat artırıyorsak, boş zamanı bir saat azaltıyor olmalıyız. Diğer eğim parametreleri de benzer bir yoruma sahiptir.

Wooldridge Veri Sorusu

Sağlık görevlilerinin ilgilendiği sorunlardan biri, hamilelik sırasında sigara içmenin bebek sağlığı üzerindeki etkilerini belirlemektir. Bebek sağlığının bir ölçüsü doğum ağırlığıdır; çok düşük doğum ağırlığı, bebeği çeşitli hastalıklara yakalanma riskine sokabilir. Doğum ağırlığını etkileyen sigara içimi dışındaki faktörlerin sigara ile ilişkili olması muhtemel olduğundan, bu faktörleri dikkate almalıyız. Örneğin, daha yüksek gelir genellikle daha iyi doğum öncesi bakıma erişimin yanı sıra anne için daha iyi beslenme ile sonuçlanır. Bunu tanımlayan bir denklem

\[ bwght=\beta_0 + \beta_1 cigs+ \beta_2 faminc + u \]

Wooldridge BWGHT data setini kullanın.

Modelin değişkenlerinin ne anlama geldiğini yazın.
Modeli tahmin etmeden, $\beta_2$ için en olası işaret nedir?
cigs ve faminc’in ilişkili olabileceğini düşünüyor musunuz? Korelasyonun pozitif mi negatif mi olabilir?
Şimdi, BWGHT’deki verileri kullanarak, faminc olan ve olmayan denklemi tahmin edin. sonuçları rapor edin. örnek boyutu ve R-kare dahil olmak üzere denklem formunda yazın. Sonuçlarınızı tartışın, faminc eklemenin cigs’in bwght üzerindeki tahmini etkisini önemli ölçüde değiştirip değiştirmediğine odaklanın.

Bu soruyu yanıtlamak için DISCRIM verilerini kullanın. Bunlar, New Jersey ve Pennsylvania’daki fast-food restoranlarındaki çeşitli ürünlerin fiyatlarına ilişkin posta kodu düzeyinde veriler ve posta kodu popülasyonunun özellikleridir. Buradaki fikir, fast-food restoranlarının siyahların daha yoğun olduğu bölgelerde daha yüksek fiyatlar talep edip etmediğini öğrenmektir. Modelimiz

\[ psoda=\beta_0 + \beta_1 prpblck+ \beta_2 income + u \]

Modelin değişkenlerinin ve prppov değişkeninin ne anlama geldiğini yazın.
ortalama prpblck ve income değerlerini standart sapmalarıyla birlikte bulun. prpblck ve income ölçü birimleri nelerdir?
Bu modeli OLS ile tahmin edin ve sonuçları, n ve R-kare dahil olmak üzere denklem biçiminde rapor edin. (Tahminleri raporlarken bilimsel gösterimi kullanmayın.) prpblck üzerindeki katsayıyı yorumlayın. Sizce ekonomik olarak büyük mü?
Basit regresyon \[ psoda=\beta_0 + \beta_1 prpblck + u \] modelini kullanarak basit regresyonu tahmin edin. Ayrımcılık etkisi income’ı kontrol ettiğiniz modele göre daha mı büyük daha mı küçük?
Gelire göre sabit fiyat esnekliğine sahip bir model daha uygun olabilir. \[ log(psoda)=\beta_0 + \beta_1 prpblck+ \beta_2 log(income) + u \]

Modelin tahmin edin ve tahminlerini raporlayın. Eğer prpblck .20 (20 yüzde puanı) artarsa, psoda’nın tahmini yüzde değişimi ne olur? (İpucu: Cevap 2.xx’dir, burada “xx”i doldurursunuz)

Şimdi prppov değişkenini kısım e’deki regresyona ekleyin. $\beta_1$’e ne olur?
log(income) ve prppov arasındaki ilişkiyi bulun. Kabaca beklediğiniz gibi mi?
Aşağıdaki ifadeyi değerlendirin: “log(income) ve prppov çok yüksek oranda ilişkili olduğundan, aynı regresyonda olmalarına gerek yoktur.”

Tek ebeveynli hanelerin öğrencilerin matematik performansı üzerindeki etkilerini incelemek için MEAPSINGLE’daki verileri kullanın. Bu veriler, 2000 yılı için güneydoğu Michigan’daki okulların bir alt kümesi içindir. Sosyo-ekonomik değişkenler, Posta kodu düzeyinde elde edilir (burada Posta kodu okulların posta adreslerine göre atanır).

math4, pctsgle, lmedinc ve free değişkenlerinin ne anlama geldiklerini yazın.
Math4’ün basit regresyonunu pctsgle üzerinde çalıştırın ve sonuçları normal biçimde rapor edin. Eğim katsayısını yorumlayın. Tek ebeveynliğin etkisi büyük mü yoksa küçük mü görünüyor?
lmedinc ve free değişkenlerini denkleme ekleyin. pctsgle üzerindeki katsayıya ne olur?
lmedinc ve free arasındaki örnek korelasyonu bulun. Beklediğiniz işaret var mı?
Imedinc ve free arasındaki önemli korelasyon varsa. Tek ebeveynliğin öğrenci performansı üzerindeki nedensel etkisini daha iyi tahmin etmek için bir tanesini regresyondan analizinden çıkarmanız gerektiği anlamına gelir mi? Açıklayın.
Görünen açıklayıcı değişkenlerin her biri için varyans enflasyon faktörlerini (VIF’ler) bulun. Hangi değişken en büyük VIF’ye sahiptir? Bu bilgi, tek ebeveynliğin matematik performansı üzerindeki nedensel etkisini incelemek için kullanacağınız modeli etkiler mi?

Wooldridge Veri Soruları Cevapları

library(wooldridge)
library(rmarkdown)
data("bwght")
paged_table(bwght)

help(bwght)

kullanacağımız değişkenlerin tanımları

faminc: 1988 family income, $1000s, 1988 aile geliri cigs: cigs smked per day while preg, hamileyken içilen günlük sigara sayısı bwght: birth weight, ounces, doğum ağırlığı, ons

Bir yandan, gelirdeki bir artış genellikle gıda tüketimini arttırır ve sigara ile aile geliri arasında pozitif bir ilişki olabilir. Öte yandan, daha fazla eğitime sahip ailelerin aile gelirleri de daha yüksektir ve daha fazla eğitim ile sigara içme arasında olumsuz bir ilişki vardır.
Sigara ve faminc arasındaki örnek korelasyonu yaklaşık -0.173’tür ve negatif bir korelasyona işaret eder.

ilkreg <- lm(bwght~ cigs,data = bwght)
ikincireg<- lm(bwght~ cigs+faminc,data = bwght)
library(stargazer)
stargazer(list(ilkreg,ikincireg),type = "text")

## 
## =====================================================================
##                                    Dependent variable:               
##                     -------------------------------------------------
##                                           bwght                      
##                               (1)                      (2)           
## ---------------------------------------------------------------------
## cigs                       -0.514***                -0.463***        
##                             (0.090)                  (0.092)         
##                                                                      
## faminc                                               0.093***        
##                                                      (0.029)         
##                                                                      
## Constant                   119.772***               116.974***       
##                             (0.572)                  (1.049)         
##                                                                      
## ---------------------------------------------------------------------
## Observations                 1,388                    1,388          
## R2                           0.023                    0.030          
## Adjusted R2                  0.022                    0.028          
## Residual Std. Error    20.129 (df = 1386)       20.063 (df = 1385)   
## F Statistic         32.235*** (df = 1; 1386) 21.274*** (df = 2; 1385)
## =====================================================================
## Note:                                     *p<0.1; **p<0.05; ***p<0.01

Denklem şeklinde rapor

\[ bwght=116.974-0.463 cigs + 0.093faminc \]

\[ bwght=119.772-0.514 cigs \]

Regresyona faminc eklendiğinde sigara içmenin etkisi biraz daha az ama aradaki fark çok büyük değil. Bunun nedeni, cigs ve faminc’in çok ilişkili olmaması ve faminc üzerindeki katsayının pratik olarak küçük olmasıdır. (Faminc değişkeni binlerle ölçülür, yani 1988’de 10.000$ daha fazla gelir, öngörülen doğum ağırlığını yalnızca .93 ons artırır.)

data("discrim")
paged_table(discrim)

help(discrim)

Gördüğünüz gibi veri setinde bazı değişkenler için bazı gözlemler NA değerine sahip. NA (not available) o gözlem için mevcut değil anlamına geliyor. Örneğin pfries (price of small fries, küçük kızartmanın fiyatı) değişkeninin beşinci gözlemi veri setinde bulunan beşinci restoranının küçük kızartma fiyatını bilmediğimiz anlamı taşıyor. Mevcut olmayan gözlemler analizlerde her zaman sorunlar çıkarır. Bu mevcut olmayan değişkenlerle nasıl başa çıkacağımızı yavaş yavaş öğrenmemiz gerekir.

Modelin değişkenlerinin ve prppov değişkeninin anlamı.

psoda: price of medium soda, 1st wave, orta sodanın fiyatı. prpblck: proportion black, zipcode, restoranın bulunduğu bölgede siyahi oranı income: median family income, zipcode, restoranın bulunduğu bölgenin medyan (ortanca) aile geliri. prppov: proportion in poverty, zipcode, restoranın bulunduğu bölgede yoksulluk oranı

mean(discrim$prpblck)

## [1] NA

sd(discrim$prpblck)

## [1] NA

mean(discrim$income)

## [1] NA

sd(discrim$income)

## [1] NA

Bildiğimiz mean ve sd fonksiyonlarını kullanarak ortama ve standart sapma değerlerini bulamadık. Çıkan NA sonucu bize bu değişkenlerin içinde bazı gözlemlerin mevcut olmadığını gösteriyor olabilir. discrim veri setinde 410 gözlem olduğundan her bir gözlemi kontrol edemiyorz ve bu değişkenlerin içinde kaç tane gözlemin mevcut olmadığını çıkaramıyoruz. R bize bu konuda is.na fonksiyonu ile yardımcı oluyor. is.na aslında sorduğumuz ingilizce bir soru ve is na? derken R’a mevcut olmayan gözlem var mı diye soru soruyoruz. R’da bize her bir gözlem için o gözlemin değeri olup olmadığını TRUE (doğru) ve FALSE (yanlış) olarak geri veriyor. Örnek vermek gerekirse konsola is.na(discrim$income) yazarsanız, R income değişkeninin her bir gözlemi için TRUE ve FALSE değerlerini yükler. TRUE o gözlemin mecut olmadığını, FALSE gözlemin olduğunu söylemektedir. Ancak tekrar aynı sorunla karşı karşıya kalırız. 410 gözlem için TRUE ve FALSE değerlerini okumamız ve FALSE olanları toplamamız kaç gözlemin var olmadığını bulmamıza yol açacaktır. sum(is.na(discrim$income)) income değişkenin içindeki var olmayan gözlem sayısını toplayacak ve bize verecektir.

sum(is.na(discrim$prpblck))

## [1] 1

sum(is.na(discrim$income))

## [1] 1

Gördüğünüz gibi hem prbblck hem income değişkenlerinin birer gözlemi boş değere sahip. Bu yüzden mean ve sd fonksiyonlarının NA gözlemlerinine sahip olduğunu söylememiz lazım.

mean(discrim$prpblck,na.rm = TRUE)

## [1] 0.1134864

sd(discrim$prpblck,na.rm = TRUE)

## [1] 0.1824165

mean(discrim$income, na.rm = TRUE)

## [1] 47053.78

sd(discrim$income, na.rm = TRUE)

## [1] 13179.29

fonksiyonun içine yazdığımız na.rm (na remove, çıkar) öevcut olmayan gözlemleri hesaplamadan çıkarmamızı söyler. prbblck değişkeninin ortalaması 0.11, standart sapması 0.18, income değişkeninin ortalaması 47053, standart sapması 13179 olacaktır.

Diyelim ki siz bütün değişkenler için kaç tane gözlemin mevcut olmadığını, kaç tane gözlemin var olduğunu, ortalamasını ve standart sapmasını görmek istiyorsunuz. Bu durumda vtable paketi size yardımcı olacaktır. Aşağıdaki komutu kullanmak için vtable paketini yüklemeniz gerektiğini unutmayın.

library(vtable)

## Loading required package: kableExtra

sumtable(discrim, summ=c('notNA(x)', 'countNA(x)', 'mean(x)','sd(x)'),out='return')

##    Variable NotNA CountNA       Mean        Sd
## 1     psoda   402       8      1.045     0.089
## 2    pfries   393      17      0.922     0.106
## 3   pentree   398      12      1.322     0.643
## 4    wagest   390      20      4.616     0.347
## 5     nmgrs   404       6       3.42     1.018
## 6     nregs   388      22      3.608     1.244
## 7   hrsopen   410       0     14.439      2.81
## 8       emp   404       6     17.622     9.423
## 9    psoda2   388      22      1.045     0.094
## 10  pfries2   382      28      0.941     0.109
## 11 pentree2   386      24      1.354      0.65
## 12  wagest2   389      21      4.996     0.253
## 13   nmgrs2   404       6      3.484      1.14
## 14   nregs2   388      22      3.608     1.244
## 15 hrsopen2   399      11     14.466     2.752
## 16     emp2   397      13     17.567     8.607
## 17  compown   410       0      0.344     0.476
## 18    chain   410       0      2.117      1.11
## 19  density   409       1   4561.803  5132.408
## 20   crmrte   409       1      0.053     0.047
## 21    state   410       0      1.193     0.395
## 22  prpblck   409       1      0.113     0.182
## 23   prppov   409       1      0.071     0.067
## 24  prpncar   409       1      0.115     0.117
## 25   hseval   409       1 147399.267 56070.468
## 26  nstores   410       0      3.139     1.809
## 27   income   409       1  47053.785 13179.286
## 28   county   410       0     13.659     8.045
## 29   lpsoda   402       8       0.04     0.085
## 30  lpfries   393      17     -0.088     0.115
## 31  lhseval   409       1     11.829     0.389
## 32  lincome   409       1      10.72     0.284
## 33 ldensity   409       1      7.959     0.996
## 34       NJ   410       0      0.807     0.395
## 35       BK   410       0      0.417     0.494
## 36      KFC   410       0      0.195     0.397
## 37       RR   410       0      0.241     0.428

discrimreg <- lm(psoda~prpblck+income, data = discrim)
summary(discrimreg)

## 
## Call:
## lm(formula = psoda ~ prpblck + income, data = discrim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.29401 -0.05242  0.00333  0.04231  0.44322 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 9.563e-01  1.899e-02  50.354  < 2e-16 ***
## prpblck     1.150e-01  2.600e-02   4.423 1.26e-05 ***
## income      1.603e-06  3.618e-07   4.430 1.22e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08611 on 398 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.06422,    Adjusted R-squared:  0.05952 
## F-statistic: 13.66 on 2 and 398 DF,  p-value: 1.835e-06

\[ psoda=0.956 + 0.115 prpblck+ 0.0000016 income + u \]

Örnek boyutu 399 gözlemdir (398 serbestlik derecesi ve 9 eksik gözlem ile gösterilir) ve ayarlanmış $R^2$ 0.595’tir. prpblck katsayısı, her şey eşit olduğunda, prpblck %10 artarsa, soda fiyatının ekonomik olarak önemli olmayan derecede yaklaşık 1,2 sent artacağını gösterir.

basitdiscrimreg <- lm(psoda~prpblck, data = discrim)
summary(basitdiscrimreg)

## 
## Call:
## lm(formula = psoda ~ prpblck, data = discrim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.30884 -0.05963  0.01135  0.03206  0.44840 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.03740    0.00519  199.87  < 2e-16 ***
## prpblck      0.06493    0.02396    2.71  0.00702 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0881 on 399 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.01808,    Adjusted R-squared:  0.01561 
## F-statistic: 7.345 on 1 and 399 DF,  p-value: 0.007015

Basit regresyon ile prpblack üzerindeki katsayının tahmini 0.065’tir. Bu, önceki tahminden daha düşüktür ve bu nedenle, gelir hariç tutulduğunda ayrımcılık etkisinin azaldığını gösterir.

logdiscrimreg <- lm(log(psoda)~prpblck+log(income), data = discrim)
summary(logdiscrimreg)

## 
## Call:
## lm(formula = log(psoda) ~ prpblck + log(income), data = discrim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.33563 -0.04695  0.00658  0.04334  0.35413 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.79377    0.17943  -4.424 1.25e-05 ***
## prpblck      0.12158    0.02575   4.722 3.24e-06 ***
## log(income)  0.07651    0.01660   4.610 5.43e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0821 on 398 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.06809,    Adjusted R-squared:  0.06341 
## F-statistic: 14.54 on 2 and 398 DF,  p-value: 8.039e-07

paste( (0.2*100)*0.122, "yüzdelik artış")

## [1] "2.44 yüzdelik artış"

“Prpblck” yüzde 20 artarsa, psoda tahmini olarak %2,44 artacaktır.

logdiscrimregprpov <- lm(log(psoda)~prpblck+log(income)+prppov, data = discrim)
summary(logdiscrimregprpov)

## 
## Call:
## lm(formula = log(psoda) ~ prpblck + log(income) + prppov, data = discrim)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.32218 -0.04648  0.00651  0.04272  0.35622 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.46333    0.29371  -4.982  9.4e-07 ***
## prpblck      0.07281    0.03068   2.373   0.0181 *  
## log(income)  0.13696    0.02676   5.119  4.8e-07 ***
## prppov       0.38036    0.13279   2.864   0.0044 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08137 on 397 degrees of freedom
##   (9 observations deleted due to missingness)
## Multiple R-squared:  0.08696,    Adjusted R-squared:  0.08006 
## F-statistic:  12.6 on 3 and 397 DF,  p-value: 6.917e-08

prppov eklemek, prpblck katsayısının 0,0738’e düşmesine neden olur.

cor(log(discrim$income), discrim$prppov, use = "complete.obs")

## [1] -0.838467

Korelasyon yaklaşık olarak -0.838’dir. Bu mantıklı, çünkü gelirdeki düşüşlerin daha yüksek yoksulluk oranlarıyla sonuçlanması beklenebilir.

Yüksek düzeyde ilişkili olmalarına rağmen, her ikisinin de dahil edilmesi mükemmel bir doğrusallık ile sonuçlanmaz ve bunun yerine, ayırt edici etkiyi izole etmeye yardımcı olan başka bir kontrol değişkeni ekleyerek modeli tamamlar.

data("meapsingle")
paged_table(meapsingle)

help(meapsingle)

math4, pctsgle, lmedinc ve free değişkenlerinin ne anlama geldikleri

math4: percent satisfactory, 4th grade math, matematik başarı yüzdesi pctsgle: percent of children not in married-couple families, evli-çift ailelerde olmayan çocukların yüzdesi lmedinc: log(medinc), medinc: zipcode median family, $ (1999), bölgenin ortanca geliri free: percent eligible, free lunch, bedava öğle yemeğine uygun görülen yüzdesi

basitreg3<- lm(math4~pctsgle, data = meapsingle)
summary(basitreg3)

## 
## Call:
## lm(formula = math4 ~ pctsgle, data = meapsingle)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -47.791  -8.310   1.600   8.092  50.317 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 96.77043    1.59680   60.60   <2e-16 ***
## pctsgle     -0.83288    0.07068  -11.78   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 12.48 on 227 degrees of freedom
## Multiple R-squared:  0.3795, Adjusted R-squared:  0.3768 
## F-statistic: 138.9 on 1 and 227 DF,  p-value: < 2.2e-16

coklureg3<- lm(math4~pctsgle+lmedinc+free, data = meapsingle)
summary(coklureg3)

## 
## Call:
## lm(formula = math4 ~ pctsgle + lmedinc + free, data = meapsingle)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -34.919  -7.195   0.931   7.313  50.152 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 51.72322   58.47814   0.884    0.377    
## pctsgle     -0.19965    0.15872  -1.258    0.210    
## lmedinc      3.56013    5.04170   0.706    0.481    
## free        -0.39642    0.07035  -5.635  5.2e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.7 on 225 degrees of freedom
## Multiple R-squared:  0.4598, Adjusted R-squared:  0.4526 
## F-statistic: 63.85 on 3 and 225 DF,  p-value: < 2.2e-16

cor(meapsingle$free,meapsingle$lmedinc)

## [1] -0.7469703

library(car)
vif(coklureg3)

##  pctsgle  lmedinc     free 
## 5.740981 4.118812 3.188079

3. Bölüm: R ile Çoklu Regresyon Analizi: (Wooldridge Örnekleri) Ders Notu

Dr. Hüseyin Utku Demir

11/24/2021

3. Bölüm: R ile Çoklu Regresyon Analizi: (Wooldridge Örnekleri) Ders Notu

İki bağımsız değişkenli model

k bağımsız değişkenli model

Regresyon Tahmin Değerleri ve Hata Payı (OLS Fitted Values and Residuals)

Çoklu regresyon tahmin değerleri ve hata payının özellikleri.

Çoklu Regresyon Modeli (Multiple Linear Regression (MLR) Model) Varsayımları

Gözlemlenemeyen Değişken Sapması (Omitted Variable Bias)

Çoklu regresyon (En küçük kareler (OLS)) tahmincilerinin varyansı

Homoskedastisite (Eş varyanslık), Çoklu eşdoğrusallık (Multicollinearity) ve Varyans enflasyon faktörü (Variance inflation factor (VIF))

Ödev

Wooldridge Bölüm Soruları

Cevaplar

Wooldridge Veri Sorusu

Wooldridge Veri Soruları Cevapları